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Abstract 

A common approach in model reduc- 
tion is balanced truncation, which is based 
on gramian matrices classifying certain at- 
tributes of states or parameters of a given 
dynamic system. Initially restricted to lin- 
ear systems, the empirical gramians not only 
extended this concept to nonlinear systems, 
but also provide a uniform computational 
method. This work introduces a unified soft- 
ware framework supplying routines for six 
types of empirical gramians. The gramian 
types will be discussed and applied in a 
model reduction framework for multiple- 
input-multiple-output (MIMO) systems. 



1 Introduction 




In a control system setting, balanced truncation is 
a well known technique for model reduction. Intro- 
duced by PQ, gramian matrices were employed to 
determine controllability and observability of linear 
systems. From these gramians a balancing transfor- 
mation can be computed, enabling the truncation, 
for example, of states that are neither controllable 
nor observable. 

With [2], empirical (controllability and observabil- 
ity) gramians were introduced, which correspond 
to the analytical gramians for linear systems, while 
extending the concept of gramians to nonlinear sys- 
tems and a generally given by: 

x{t) = f{x,u,p) 
y(t) = g(x,u,p), 

with the system function / and output function 
g of states x, input u and parameters p. In the 
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special case of an unparameterized linear system 
/ = Ax(t) + Bu(t) and g = Cx(t). These empirical 
gramians are computed by averaging simulations or 
experimental data with perturbations in inputs and 
initial states. 

The emgr framework presented here, encompasses 
six empirical gramians, namely the controllability, 
observability, cross, sensitivity, identifiability and 
joint gramian. To adapt the computation of empir- 
ical gramians to the operating setting of the sys- 
tem, the initial state and the input are the main 
parameters which are perturbed by rotations and 
scaling. The sets of rotations provided are {1} (unit 
matrix) and {—1, 1} (negative unit matrix and unit 
matrix). Though these are very basic sets, and thus 
might not reflect all dynamics, especially with in- 
terrelated states and parameters, they allow a very 
efficient gramian assembly. Scales may be freely 
chosen. The subdivision of the scales may be lin- 
ear, logarithmic or geometric. Finally, there are 
several options to average against: the arithmetic 
average ([2]), a steady state (|3|), and additionally 
the principal components of the simulations or data 
via a principal component analysis (PCA). 

2 Empirical Gramians 

Concerned with the reduction of states, the con- 
trollability, observability and cross gramian are pre- 
sented next; followed by the sensitivity, identifiabil- 
ity and joint gramian, which are used for parameter 
and combined reduction. For the purpose of defin- 
ing the gramians, a linear time-invariant control 
system is assumed: 

±{t) = Ax(t) + Bu(t) 
y(t) = Cx(t), 

with the states x e R™, input u 6 R rn , system ma- 
trix A £ R nx ™, input matrix B e R nxm and output 
matrix C E R° . 



Necessary perturbations are given by six sets, of 
which with index u define the input perturbation, 
while sets those with index x define the initial state 
perturbation: 





= {ei G [ 


? m ;|| ei || =l;e iej 


= 0;i = 


1, . . . ,m 


E x 


= {fi e I 


R n ;ll/i|| = i;/i/j 


= 0;i = 


l,...,n} 


Ru 


= {R G 


R mxm ; = 1 


; i = 1, . 


..,t} 


Rx 


= {S> G 


R nx "; SfS'i = l; t 


= 1,... 


,*} 


K u 


= {h G 1 


R; fc 4 > 0; i = 1, . . 


•,«} 




K x 


= {Zi G F 


l;l i >0;i=l ) ... 


,[}■ 





These sets should correspond to the ranges in in- 
puts and initial states the system is operating in. 

2.1 Controllability Gramian 

Controllability is a quantification of how well a 
state can be driven by input. Analytically, the con- 
trollability gramian is given by the smallest semi- 
positive definite solution of the Lyapunov equation: 
AW C + W C A T = -BB T . If the underlying system 
is asymptotically stable the controllability gramian 
can also be defined using the linear input-to-state 
map: 

POO 

W c = J e AT BB T e ATT dr. 

Following |3] , the empirical controllability gramian 
is defined by: 

Definition 

For sets E u , R u , K u , input v(t) and input during 
the steady state x, u, the empirical controllabil- 
ity gramian is given by: 



r t m 



^^EEE^/ *^ 

n h=l i= l j=l K h Jo 
^ hlj (t) = (x hij (t) - x)(x hij (t) - x)* G R nXTl . 

With x h% i being the states for the input configura- 
tion u hl i{t) = khRiejv(t) + u. 

Originally, in [2J, v(t) was restricted to 5(t), but 
extended in [3J to arbitrary input under the name 
empirical covariance matrix, x can be arithmetic 
average or steady state. Restricting E u to { — 1,1} 
simplifies the input perturbation to: 

u Mj (t) = -Vk h e ]V (t) + u. 



2.2 Observability Gramian 

Observability quantifies how well a change in a 
state is reflected by the output. The analytical ob- 
servability gramian is given by the smallest semi- 
positive definite solution of the Lyapunov equation: 
AWo + W A T = -C T C. Given an asymptotically 
stable underlying system the observability gramian 
can also be defined using the state-to-output map: 



W = I e A T C T Ce AT dT. 



The empirical observability gramian is defined as 
described in [2] and [3J: 

Definition 

For sets E x , R x , K x and output during the steady 
state x, y, the empirical observability gramian 

is given by: 
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K r b = (y qrs (t)-yr(y qrs (t)-y)£ 



(i) 



With y qrs being the systems output for the initial 
state configuration x^ rs = l q S r f s + x. 

y can be arithmetic average or steady state output . 
Restricting E x to {—1,1} simplifies (*) to: 



Wr. 
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and the initial state perturbation to: 
-l%f s + x. 



(IIS 



2.3 Cross Gramian 

The cross gramian ([?]) makes a combined state- 
ment about the controllability and observability, 
given the system has the same number of inputs 
and outputs. If the system is also symmetric, mean- 
ing the system transfer function is symmetric, then 
this gramians' eigenvalues equal the Hankel singu- 
lar values. It is originally computed as the small- 
est semi-positive definite solution of the Sylvester 
equation: AW X + W X A T = -BC. The cross 
gramian can also be defined using the input-to-state 



and state-to-output maps, if the underlying system 
is asymptotically stable: 



W x = 



z Ar BCe A T dr 



The empirical cross gramian has been introduced in 
[5] for SISO systems and was extended to MIMO 
systems in [5]. 

Definition 

For sets E u , E x , R u , R x , K u , K x , input u during 
steady state x with output y, the empirical cross 
gramian is given by: 

h=l i=l 9=1 r=l H u 

9** = fJSjAx'^ (t)ejRfAy^ (t) (**) 
Ax hij {t) = (x hij (t) - x) 

Where x hl i and y qrs being states and output for 
the input u h ^(t) = khRiejv(t) + u and initial state 



qrs 



IqSrfs + x respectively. 



x and y can be for example arithmetic average or 
steady state (output). Again, restricting E u and 
E x to {—1, 1} simplifies (**) to: 



hiqr 



(-iy + \fjAx h ^(t)ejAy^(t), 



as well as input and initial state perturbation to: 



u hij (t) = -Vkhejvit) +u 

Xq 1 lq f s ~l~ X . 

2.4 Sensitivity Gramian 



The sensitivity gramian allows controllability based 
parameter reduction and identification. It is based 
on [6] and aimed for models that can be partitioned 
as follows: 

p 

x = f(x,u,p) = f(x,u) + ^2f{x,p k ). 

fe=i 

The parameters p G R p are handled here as ad- 
ditional inputs. All summands of the partitioned 
system are treated as independent subsystems and 
thus a controllability gramian for each subsystem 



can be computed. Each parameters controllability 
is encoded in the sum of singular values of the as- 
sociated (sub-) controllability gramian Wc,k- The 
sensitivity gramian is now given by the diagonal 
matrix, with each diagonal element being the trace 
of a sub controllability gramian: 

fc=i 



The controllability and thus identifiability of each 
parameter is then given by the corresponding di- 
agonal entry of the sensitivity gramian Ws- For 
partionable linear systems, the sum of all subsys- 
tems controllability gramians Wc.k and the param- 
eter free subsystems gramian Wc,o equals the usual 
controllability gramian (0): 

> 

The sensitivity gramian can be applied to non- 
partionable models with reduced accuracy. 

2.5 Identifiability Gramian 

The identifiability gramian enables observability 
based parameter reduction. As described in [TJ, 
the dynamic systems states are augmented with as 
many states as parameters that are constant over 
time, and have the initial value of the (prior) pa- 
rameter value. 

'•A ( f(x,u,p) 

y = g(x,u,p) 

The observability gramian of this augmented sys- 
tem holds the observability information of states 
and parameters. To extract the parameter specific 
observability, the Schur-complement can be applied 
to the overall gramian: 
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^Wi = W P - WqWo~ x W q S. 
2.6 Joint Gramian 

Based on the identifiability gramian procedure, the 
cross gramian can be employed for a concurrent 



state and parameter reduction (see |5]). Not only 
augmenting the states with constant parameter 
states as before, but also adding as many inputs 
v and outputs as parameters, acting via identity 
on the augmented states and augmented outputs 



the empirical observability gramian assembly from 
equation Q] can be computationally simplified to: 
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to preserve symmetry. The cross gramian of this 
special augmented system, similar to the identifia- 
bility gramian, holds the cross gramian of the origi- 
nal system as well as a cross identifiability gramian 
Wj which can be extracted with Schur-complement 
from the joint gramian: 



Computation of empirical gramians using emgr 
is very portable, since only basic vector and ma- 
trix operations are required. Necessary integra- 
tions, meaning simulations for given inputs or ini- 
tial states, are accomplished either by the first or- 
der Euler's method or the second order Adams- 
Bashfort method. The empirical gramian frame- 
work emgr as well as the following experiments 
are released under an open source license, are com- 
patible with Octave and Matlab, and can be found 







at: http://gramian.de . 
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Wj 

=^Wj = W P - WqWx^Wq. 

3 Implementation 



The emgr software framework presented here pro- 
vides a uniform interface to compute all six empir- 
ical gramians and is given by: 

W = emgr(f , g, q, p, t, w, v cfg , u, u s , x s , u m , x m ); 

with f and g being handles to the system and the 
output function both with requiring the signature 
f(x,u, p) and g(x,u, p). q is a vector defining the 
systems number of inputs, states and outputs; p 
holds the parameters and u provides the input to 
f and g; t is a three component vector containing 
start time, time step and stop time; w is a charac- 
ter setting the gramian type. The ten component 
vector v cfg configures the available options, includ- 
ing averaging types, symmetry enforcement, input 
and state scale subdivisions and perturbation rota- 
tions. x s and u s set a steady state and the associ- 
ated input. Lastly, u m and x m define the scales of 
the perturbation. 

The parameter reducing empirical gramians (sen- 
sitivity, identifiability and joint) are an encapsula- 
tion of the state reducing empirical gramians (con- 
trollability, observability and cross). Computation 
of the latter is extensively vectorized, exploiting 
the gramian matrix assembly format. In example 



4 Numerical Experiments 

To demonstrate the various empirical gramians, 
computed by the emgr framework, a symmetric 
nonlinear MIMO system with a hundred states, ten 
inputs and ten outputs is employed. The system 
matrix is generated randomly with ensured sym- 
metry; the input matrix B is also a random matrix 
and the output matrix is given by C = B T . Fur- 
thermore a random parametrized source term E p is 
added. Input is applied through a delta impulse. 



f(x,u,p) 
#0) 



y 



Ax 3 
Cx 



Bu + E n 



For the unreduced system output see figure [T] 
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Figure 1: Original system output. 

First, a state reduction, using the empirical con- 
trollability gramian and the empirical observability 
gramian, through balanced truncation is performed 
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the input is applied, yet yield overall qualitatively 
similar results. The empirical sensitivity gramian 
can be applied if the underlying system can be par- 
titioned such that f(x,u,p) = f(x,u) + f(x,p). To 
be able to use it in this setting, the source term is 
parametrized with as many parameters as states, 
which is reduced to the number of outputs in figure 

SI 




90 100 




Figure 2: Reduced system output by balanced trun- 
cation using the empirical controllability gramian 
and the empirical observability gramian WC = 
emgr(f ,g, [10, 100, 10] ,p, [0,0 .01 , 1] , 'c' ) ; W0 
= emgr(f ,g, [10, 100, 10] ,p, [0,0. 01 , 1] , 'o ') ; 



in figure [3J reducing the number of states to the 
number of outputs. 

Next, a state reduction by truncation employing 
the empirical cross gramian is demonstrated in fig- 
ure [3j again reducing the number of states to ten. 
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Figure 4: System with reduced source 
term; reduction by truncation using 
the empirical sensitivity gramian WS = 
emgr(f ,g, [10 , 100 , 10] ,p, [0, . 01 , 1] , ' s ' ) ; 

Since the parameters of the source term are re- 
duced, the cumulative effects in the original sys- 
tem are the origin of the increasing error over time. 
Next, the parametrized source term is reduced by 
the empirical identifiability gramian in figure [3] 
Finally, in figure the same system undergoes a 
combined state and parameter reduction using the 
empirical joint gramian. The error here is due in 
part due to the combined reduction of first param- 
eters and second states as this affects the states 
twice. 

5 Future Work 



Figure 3: Reduced system output by trunca- 
tion using the empirical cross gramian WX = 
emgr(f ,g, [10,100,10] ,p, [0,0.01,1] , 'x') ; 

Both approaches produce a peak in the error when 



The emgr framework already allows a wide range 
of computations of empirical gramians for state or 
parameter reduction. 

Further work will enhance the flexibility, while 
keeping the interface as simple as possible. Fol- 
lowing [TJ, allowing factorial designs will greatly 
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